#if !defined(IMGLIB_GRAD_H)
#define IMGLIB_GRAD_H

#include "typelib/axis.h"
#include "imglib/image.h"

namespace img {	

class axis: public types::axis<float> {
	typedef types::axis<float> superclass;
public:
	axis(float x = 0, float y = 0)
		: superclass(x, y) { }
	axis(point<float> p)
		: superclass(p) { }
};

} //namespace img

//same as in types

template <>
inline void inverse(img::axis & a) 
	{ inverse(static_cast<types::axis<float> &>(a)); }

template <>
inline double scalar_product(const img::axis & a0, const img::axis & a1) {
	const types::axis<float> & p0(a0), & p1(a1);
	return scalar_product(p0, p1);
}

//linear & bilinear differ from types in normalization

template <class F>
inline img::axis linear(const img::axis & a0, const img::axis & a1, F a) { 
	const types::axis<float> &t0(a0), &t1(a1);
	return img::axis(linear(t0, t1, a)).normalize();
}

template <class F>
inline img::axis bilinear(const img::axis & a00, const img::axis & a01, const img::axis & a10, const img::axis & a11, F a0, F a1) {
	const types::axis<float> &t00(a00), &t01(a01), &t10(a10), &t11(a11);
	return img::axis(bilinear(t00, t01, t10, t11, a0, a1)).normalize();
}

namespace img {	

template <class T> void gradient(const T & dx, const T & dy, point<float> * dir, float * equiv = 0, float * min = 0, float * max = 0) {

	double A = scalar_product(dx, dx);
	double B = scalar_product(dx, dy);
	double C = scalar_product(dy, dy);

	const float eps = std::numeric_limits<float>::epsilon();

	double a = A - C;
	double b = 2 * B;
	double c = A + C;
	double d = _hypot(a, b);
	double e = c - d;

	double sqrt2 = sqrt(0.5);

	if (fabs(e) <= eps)
		e = 0;

	if (c > eps) {
		if (d > eps) {
			if (dir != 0) {
				dir->x = float(sqrt2 * sqrt(1 + a / d));
				dir->y = float(sqrt2 * sqrt(1 - a / d));
				if (b < 0)
					inverse(dir->x);
			}
			if (equiv != 0)
				* equiv = float(sqrt(e / (c + d)));
			if (min != 0)
				* min = float(sqrt2 * sqrt(c - d));
			if (max != 0)
				* max = float(sqrt2 * sqrt(c + d));
		} else {
			if (dir != 0)
				* dir = point<float>();
			if (equiv != 0)
				* equiv = 1;
			if (min != 0)
				* min = float(sqrt2 * sqrt(c));
			if (max != 0)
				* max = float(sqrt2 * sqrt(c));
		}
	} else {
		if (dir != 0)
			* dir = point<float>();
		if (equiv != 0)
			* equiv = 0;
		if (min != 0)
			* min = 0;
		if (max != 0)
			* max = 0;
	}
}

} // img


#endif // IMGLIB_GRAD_H
